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We examine the interactions of a black hole with a massless scalar field using a coordinate system which extends ingoing 
Eddington-Finkelstein coordinates to dynamic spherically symmetric-spacetimes. We avoid problems with the singularity by 
excising the region of the black hole interior to the apparent horizon. We use a second-order finite diff'erence scheme to solve 
the equations. The resulting program is stable and convergent and will run forever without problems. We are able to observe 
quasi-normal ringing and power-law tails as well an interesting nonlinear feature. 



The Seidel and Suen paper is also notable for the intro- 
duction of a general technique, termed causal differenc- 
ing, which ensures that, independently of the details of 
the coordinate system adopted, the difference scheme's 
numerical domain of dependence contains the physical 
domain of dependence (i.e. is causal). Alcubierre and 
Schutz have used a similar but somewhat more general 
technique which they call causal reconnection to treat the 
wave equation on an arbitrarily moving grid Scheel 
et. al. have recently used black hole excising in a study 
of gravitational collapse in Brans-Dicke gravity 0J^. 

In this paper, we again use the black hole excising tech- 
nique to examine the interactions of a black hole with 
a massless, minimally coupled scalar field in spherical 
symmetry. Unlike previous work, we use a null based 
slicing to get a coordinate system (ingoing Eddington- 
Finkelstein) which fits naturally with black hole exci- 
sion (for a Schwarzschild black hole, all variables are 
static and non-singular (appropriately smooth) every- 
where on the solution domain, and perturbations about 
Schwarzschild give perturbations of this behavior). Fur- 
ther, we introduce a modification to these coordinates 
which allows for easy tracking of the apparent horizon. 
We are able to get a second-order, convergent evolu- 
tion scheme which will stably evolve forever, and which 
shows the expected effects of quasi-normal ringing and 
the power law decay of the scalar field. We also see some 
interesting coordinate and nonlinear effects. 

The plan of the remainder of the paper is as fol- 
lows. In Se c. O , and following some unpublished previ- 
ous work |P,|l0|, we define the minimally-modified ingo- 
ing Eddington-Finkelstein (MMIEF) coordinate system 
and derive equations of motion for the gravitational and 
scalar fields which are specialized to this coordinate sys- 
tem. This derivation is based on detailed calculations of 
the Einstein-Klein- Gordon equations of motion in a gen- 
eral spherically symmetric coordinate system [ pT| which 



I. INTRODUCTION 

A longstanding goal of relativists has been the long- 
time numerical evolution of a black hole spacetime. Such 
an evolution is difficult because of the physical singular- 
ity the spacetime contains. Despite their early promise, 
it has not sufficed to use slicings (choices of time coor- 
dinate) which avoid the physical singularity. Invariably, 
such coordinate systems develop coordinate singularities 
in the vicinity of the event horizon and computation- 
ally, physical and coordinate singularities are essentially 
equally pathological. Many years ago, Unruh suggested 
that it might help to consider evolution of only the ex- 
terior of a black hole. In fact, he argued, since the black 
hole interior is, by definition, causally disconnected from 
the rest of the universe, evolution of events within the 
horizon constituted wasted computational effort. In the 
necessary absence of exact information concerning the 
spatial location of the event horizon at any instant dur- 
ing the evolution of given black hole initial data, Un- 
ruh proposed that the apparent horizon be tracked and 
used as an approximation to the true horizon. Thornburg 
developed these ideas, first in the context of generating 
initial data for n black holes (each of which could have 
arbitrary momentum and spin) and then in a pro- 
gram of research for solving the vacuum axisymmetric 
Einstein equations [^^. However, the first clearly suc- 
cessful application of the black-hole excising technique 
in a dynamical situation came with the work of Seidel 
and Suen who studied spherical evolution of the vac- 
uum (Schwarzschild) as well as a self-gravitating massless 
scalar field Extensions of this work are discussed in 
Anninos et. al. §. 
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have been summarized in Appendix A. In Sec. [11, we ex- 
amine the issue of regularity at the origin, r = 0, with the 
principal result that use of the coordinate system must 
be restricted to cases where matter never reaches r — 0. 
We follow with a detailed description of our finite differ- 
ence approximations in Sec. IV and a discussion of initial 
data in Sec. 0. We discuss the convergence and stabil- 
ity of our difference solutions in Sec. VI — the evidence 
presented there suggests that the scheme can be used 
to carry out arbitrarily long evolutions. Various physi- 
cal and coordinate effects which hav e em erged from our 
studies are discussed in sections Sec. |VI]| -Sec. XI and we 
end with some concluding remarks. 

We note that the programs used to generate the re- 
sults described below were written in RNPL (Rapid Nu- 
merical Prototyping Language), a language designed by 
the authors This language allows for the easy and 
compact expression of time-dependent systems of partial 
differential equations and the rapid trial of various finite- 
difference techniques for solving them. 



II. MINIMALLY-MODIFIED INGOING 
EDDINGTON-FINKELSTEIN COORDINATES 
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FIG. 1. The ingoing Eddington-Finkelstein slices in 
Kruskal-Szekeres coordinates. The dotted lines are constant 
t, while the dashed lines are constant r. The dark curves are 
the singularity. The diagonal lines are the horizon (r — 2M). 

We begin by recalling that the ingoing Eddington- 
Finkelstein (lEF) coordinate system (see, for example. 
Chapter 30 of is defined for static, spherically- 

symmetric spacetimes (i.e. Schwarzschild) and combines 
an areal (measures proper surface area) radial coordi- 
nate r, with a time coordinate t, which is chosen so 
that the ingoing tangent combination dt — dr is null. 
Fig. 1^ shows some slices of constant lEF time plotted 



in Kruskal-Szekeres coordinates. Notice that all of the 
slices penetrate the horizon and meet the singularity. 

In generalizing lEF, we first consider the following 
general "3-1-1" form for a time-dependent spherically- 
symmetric metric: 

ds^ = {-a^ + a^(3^) dt^ -t- 2a^(3dtdr + a^dr^ + r^b^dn'^. 

(2.1) 

Here, a,b,a, and (3 are functions of r and t, and dO? is 
the metric on the unit-sphere. We fix the spatial degree 
of coordinate freedom by introducing a "shifted" areal 
coordinate s, defined by s = r + f{t) for some as yet 
undetermined function /. With this choice, the metric 
becomes 

ds^ = (-a^ + a^p^) dt^ + 2a^pdtdr + a^dr^ + s^dn^. 

(2.2) 



Comparison with the general form (A4) yields the iden- 
tifications s = rb and 5=1 + f /r. From the general 



evolution equation for b ( All ), we have 
{rb) = -arbK% + (3 (rb)' . 



(2.3) 



Here and subsequently, overdots and primes denote par- 
tial differentiation with respect to t and r respectively. 
Solving this last equation for the shift vector component 
(or simply "the shift") /3, and noting that {rb) = /, 
(rb)' = 1, we obtain the condition on the shift which 
must be dynamically enforced to keep the metric func- 
tion s areal: 



(2.4) 



We now fix the time slicing by demanding that the ingo- 
ing combination of timelike and radial tangent vectors — 
dt — dr — be null, exactly as in the case of the original 
lEF coordinates. This gives a condition on the metric, 
namely gtt — ^gtr + Qrr = 0. Using (p. 2]) this implies 



a = ±a{l- P) . 



(2.5) 



We choose the sign so a is positive for /? < 1, that is 
a = a (1 — /3). Using this and ( ^ ) we get 



/ + saK% 



a{l-f) 



1 + saK% 1 -I- saK<^g ' 

Hence the metric takes the form 

ds^ = (2/3 - 1) dt'^ + 2a^(3dtdr + a^dr^ + s^dn^ 



(2.6) 



(2.7) 



Factoring the first three terms yields 

ds^ = ((2/3 -l)dt + dr) {dt + dr) + s^dn^, (2., 
which shows that the characteristic speeds are 
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-1,1-2/3. 



(2.9) 



We now specialize the g enera l spherically-symmetric 
equations of motion (A1C)-(A18) given in Appendix A to 
the MMIEF coordinate system. The full set of constraint 
and evolution equations for the geometric variables is: 



n^) = 



K% - K\ 47r$n 



a = -a'^ {I - 13) K\ + {a(3)' 



K% = pK^'g + a (1 - /3) K'g {K\ + 2K^ 
V a I as 



(2.10) 



(2.11) 



(2.12) 



(2.13) 



I3K^', + a (1 - /3) K% {K\ + 2K%) 



(3-1 
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/3'a' 






a 



— h 87r$^ 

a / sa 



(2.14) 



while the wave equation for the massless scalar field be- 
comes the first-order (in time) system: 

6 = (/3$ -h (1 - /?) n)' (2.15) 



U=^[s^f3n+{l-(3)<^>)]' -^U. (2.16) 

We note that we have used the slicing condition, a = 
a(l — (3), to eliminate the lapse function from the above 
set. In addition, as is always the case in general rela- 
tivistic dynamics, we have more constraint and evolution 
equations governing the geometric variables than vari- 
ables themselves. We adopt the often-used (and some- 
what ad hoc) approach of discretizing some sufficient sub- 
set of the equations with the expectation that, provided 
the difference scheme converges, residuals of discretized 
forms of the remaining equations will be of the same or- 
der in the mesh spacing as the truncation error of the 
difference scheme itself. In fact, we have made consider- 
able use of our freedom to construct schemes based on 
various combinations of constraint and evolution equa- 
tions in the development of the stable second-order meth- 
ods described in Sec It is entirely possible that we 
could construct a stable scheme without explicit use of 
the constraints but our attempts to do so with the type 
of differencing described in Sec IV were not successful. 



Since the function f{t) (recall s{r,t) = r + f{t)) is still 
unspecified, we require one more equation to determine 
the time-evolution of our model system. We derive an 
evolution equation for / by demanding that a certain 
radial coordinate, r — r/i, be tied to the location of a 
marginally trapped surface. In general we will want to 
track the outermost marginally trapped surface, so in the 
following we will assume that the surface we are tracking 
is, in fact, an apparent horizon. We recall that if S^^ is an 
outward-pointing, space-like unit normal to a marginally 
trapped surface, then it satisfies p3[ 



0, 



(2.17) 



where Di is the covariant derivative compatible with the 
3-metric and K = K^i. In spherical symmetry, and with 
the functional forms of the 3-metric and extrinsic curva- 
ture given in Appendix A, this equation reduces to 



{rh)' ~ arbK\ = 0. 
In the MMIEF system this is simply 
asK\ = 1. 



(2.18) 



(2.19) 



To keep the apparent horizon at fixed r, we demand that 
[asK^ e)\r^ — 0, where is the initial position of the 
apparent horizon. Solving this equation for / and using 
the constraint and evolution equations to eliminate time 
and space derivatives of the geometric variables, we find: 



/ 



47rs2 ($ + n) 



(2.20) 



We are most interested in the situation where the ini- 
tial data for our spacetime describes a black hole of mass 
M which is well separated from any scalar field. In this 
case we have Vh = 2M and, provided that ( ^.20 ) is satis- 
fied, the apparent horizon will remain a,t r ~ rh — 2M . 
The area of the apparent horizon, however, is given by 
^TTs\ = 4:TT{rh + /)^, and hence will increase as matter 
falls into the black hole, in accord with physical expec- 
tations. 

It is useful to write down the Schwarzschild solution in 
lEF coordinates. The metric is usually written as 

ds^ = -(^^- -y-^ dV"^ + "^dVdr + r^dn^, (2.21) 

where V is a null coordinate. Defining a time-like coor- 
dinate t = V — r, the metric becomes: 



ds' = - 1 - 



1-1- 



2M 

r 
2M 

r 



dt-^ H dtdr 

r 



dr^ + r^dn^. 



(2.22) 



Comparison of this last result with the general form ( A4) 
yields the following expressions for the various metric 
components: 
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r + 2M 



2M 
r + 2M 



r + 2M 



(2.23) 



where we use overbars on quantities to stress that the 
expressions are vahd only for the vacuum case. Using the 



above and equations (AlO) and (All), we can compute 



the extrinsic curvature components: 

K% = 1E(L±^ ^ -2Mir + M)^ ^^^^^^ 

(r(r + 2M))2 {r{r + 2M))^ 

Since we are working in spherical symmetry, we can 
meaningfully define the mass (or mass aspect) function 
rn(r, t) which, at least in a vacuum region, provides an 
invariant measure of the gravitational mass contained 
within radius r at time t. Moreover, even when mat- 
ter is present, m{r,t) and m'{r,t) are useful diagnostic 
quantities in our calculations. The mass in MMIEF co- 
ordinates can be computed from the surface area using 
the general expression: 

m (r, t) - is (l - {WnA)-^ A-'^^,^) , (2.25) 

where A — Att s'^. Substituting this expression for the 
area into ( ^.25| ) and differentiating, we get 



m (r, t) 



1^ / l-{saK%y 



(2.26) 



By making use of the evolution and constraint equations, 
we can write this mass as an integral over the mass- 
density. In this form we have 



Att I s- 



2a2 



sK\^ 1 dr. 
a 



(2.27) 



where as before, Vh is the location of the apparent hori- 
zon and Sh = s{rh). We note here that as used above, 
M represents the mass of the black hole, that is Sh/2. 
We also define M^a as the total mass in the spacetime, 
namely 



= m (cxd) 

+ 477 



M 



sK\—\dr. (2.28) 



$2 + n2 

2a2 



Given that our finite difference grid does not extend to 
infinity, we will approximate M^o by taking the upper 
limit of the above integral to be the outer boundary of 
our computational domain. 



III. REGULARITY AT THE ORIGIN 

In cases where a black hole is not initially present or 
there is insufficient mass in the scalar field to form a 



black hole via collapse, the infalling matter will encounter 
the coordinate origin, r = 0. As is generically the case 
when using spherical coordinates, the various geometric 
and matter variables must satisfy regularity conditions 
as r ^ in order that the origin remain a regular point 
in the spacetime (see for an extensive discussion of 
regularity conditions). 

Since is a scalar and a,K^g, and K^r are components 
of rank-two tensors, we assume that they are even in r 
near the origin. This means that their spatial derivatives 
must vanish at r = 0: 



0. 



(3.1) 



Elementary flatness near the origin in MMIEF coor- 
dinates dictates that a{0, t) — 1. An examination of 
the momentum constraint ( 2.11 ) shows that K^g{Q,t) = 
K^riO,t). We can find further conditions by examining 
the pote ntially divergent terms of the evolution equation 
for K^g ( ^.13 ). These terms, which are those with pow- 
ers of r in the denominator, can be collected and written 
as 



(1 - /?) (a2 - l) + r/3' 



(3.2) 



Clearly, a s r — > both the numerator and the denom- 
inator of ( ^.2[ ) go to zero, so we use L'Hopital's rule to 
compute the limit. The derivative of the numerator is 
2 (1 - /?) aa' - /?' (a2 - l) + + r/3". As r this goes 
to (3'. The derivative of the denominator is r^a' + 2ra. 
Clearly this goes to zero as r goes to zero. Thus, we must 
have P{0,tY = 0. 

Since we still have an indeterminate form, we ap- 
ply L'Hopital's rule again. The second derivative of 

the numerator is (1 - (3) (2aa" + 2 [a'f^ -13" (a^ - l) - 

2af3'a' + 2/3" -I- r/3"'. As r goes to zero, this goes to 
2 (1 — /3) a" -I- 2(3" . The second derivative of the denom- 
inator is r^a" + 4ra' -|- 2a. The limit of this is 2. Thus, 
we have 



. (1-/3) (a2 - 1) + r/3' 



lim 



{I - 13) a" + 13" . (3.3) 



Let us now consider the behavior of the spatial deriva- 
tives of /3 as r ^ 0. When there is no black hole present, 
the shift is given by 



raK" 



1 -I- raK^g ' 



(3.4) 



From this we can easily see that (3 is zero at the origin. 
Now 



_ raii'% + ra'K\ + aK'^ 



{l + raK\ 



Thus, 



lim B' 



(3.5) 



(3.6) 
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Since we have argued that (3'{0,t) = and K^g{0,t) 
K^r{0, t), we have 



(3.7) 



K''eiO,t)=K%{0,t) = 0. 
Now the second derivative of /3 is 

raK^'' + 2aK^' + 2ra'K^' + 2a' K\ + ra" K% 



{l + raK%y 



2 ((raX%)' 
{l + raK%) 



(3.8) 



As r 0, this expression vanishes. Thus, P" (0, t) = 0. 

Since K'^g{Q,t) = 0, we must ha ve k \{Q,t) = and 
consequently the right hand side of ( 2.13| ) must vanish as 
r ^ 0. This wiU happen only if the limit of ( |3.2[ ) is zero, 
which, given the results deduced above, can only happen 
if a(0,t)" = 0. 

The evolution of the scalar field is accomplished 
through tw o aux iliary variables, (f> and 11. These are 
defined b y ( | A14| ). The condition on $ has already been 
stated in ( |3.l| ) 



*(o,t) = o. 



(3.9) 



and it is easy to show that 11 must be even in r near the 
origin. Thus we have 



n'(0,i) = 0. 



(3.10) 



TABLE I. Two-Level Finite Difference Operators 



Operator 


Definition 




Expansion 




(-3/r + 4/;Vi - /r+2) /2Ar 


a./ 


" + 0(Ar2) 


A^/r 


(3/r-4/r_i + /r_2)/2Ar 


a./ 


'^ + 0(Ar2) 


A./r 


(/r+i-/r_i)/2Ar 


a./ 


Vo(Ar^) 
' + O (At^) 


At/r 


(/r+' - /r) /At 




A?/r 


- /r) /At+ 


a/r+'+o(At2) 












4(/r_i + /;;i)]/i6At 








(/r+' + /r) /2 


ft 


^ + O (At^) 




(/r + /;ii)/2 


1 + 0(Ar2) 


A^vr 


MtA//r 


O (A> + At^) 




/^tAi^/r 








o 


(Ar2 + At') 




MtAr/r 








o 


(Ar^ + At^) 


A?/r 








/r+i-/r)/2Ar 


O (Ar^ + At^ 








+ArAt) 



Unfortunately, these regularity conditio ns ar e incon- 
sistent with the Hamiltonian constraint ( 2.10 ). If we 
solve ( ^.10| ) for a', take a radial derivative of the resulting 
expression, then take the r ^ limit, we find: 



a"(0,i) oc n(0,i)2 ^ 0. 



(3.11) 



Thus, we are lead to the conclusion that the MMIEF co- 
ordinate system will admit no non-singular curvature at 
the origin. The only consistent solutions near the origin 
describe fiat space or a black hole. Thus, MMIEF is a 
"good" coordinate system to use only when a black hole 
already exists in the spacetime. In a spacetime without 
a black hole, the equations will remain consistent as long 
as no matter encounters the origin. This will be the case 
if the scalar field is outgoing or if it collapses to form a 
black hole before it encounters the origin. For a collapse 
problem, we could start with another coordinate system 
and change to MMIEF coordinates if an apparent horizon 
forms. If no such horizon forms, there is really no need for 
the special horizon tracking properties of MMIEF coor- 
dinates anyway. Though we have not implemented such 
a scheme, this approach was used successfully in @J|,0,||. 



IV. FINITE DIFFERENCE EQUATIONS 



We use equations (^T2|) , ( ^T3|) , ( |2T5| ), (|2T6|) , and 
( ^.20D to evolve a, K%, g^^and /; equation ( |2.1l| ) 
to find K^r\ and equation ( P^ ) to find j3. We solve 
these equations using finite difference techniques on a 
uniform mesh with spacings Ar and At = AAr, where 
the Courant factor A, is held fixed when we change the 
basic scale of discretization. 

Table || shows the operators we use in the discretiza- 
tions. Note that while the derivative operators take a 
lower precedence than the arithmetic operators, that is 
^rfi^'^ = (/i+i^ — /i-i^) /Ar, the time averaging oper- 
ator takes a higher precedence, that is fitfr^ = {fJ-tfr)'^ 
andMi(a«6f)^=Mia>,6f. 

In the interior, the finite difference equations are: 



(a^ (1 -/?));' + A^ (a/3): 



(4.1) 



Aj K = ml3^ A^K e, + fJ-t { — — [a-- 

A a an 

+ f^t (a (1 - P) K\ {2K\ + K\)yi , (4.2) 



s « / , 



A^ (/3<i> + (i-/3)n):, 



(4.4) 
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rA^ {s^ (/3n+(l 



and 



(4.5) 



(4.6) 



(4.7) 



(4.8) 



These equations are applied everywhere in the interior 
except at the two points next to the boundary points. 
At these points, we use the same equations except the 
dissipative time derivatives (A^) are replaced by regular 
time derivatives (At), since the value at i + 2 or i — 2 is not 
available at these locations. It is interesting to note that 
all of the spatial derivatives ar e a ngled (A^) except for 
the derivative of (3 in equation (12) and the derivative of 
B in equation (4.3). Switching any of these derivatives 
from angled to non-angled or from non-angled to angled 
results in an instability. We don't have an explanation 
for why this particular combination of derivatives works, 
but the use of angled derivatives here was motivated by 
their successful application in previous work by Choptuik 
(see and ) . This is a good example of a technique 
that was perfected by "experimental" numerical analysis. 

As discussed previously, the inner boundary of our 
computational domain is fixed to the apparent horizon, 
r — Th- In accord with the causal properties of the ap- 
parent horizon, function values on the horizon can be ad- 
vanced without use of any information defined a,t r < r^- 
Therefore, dX r = Vh we use the same equations employed 
in the interior, except that we replace centered differences 
with forward differences. For example, equation (4.1) be- 
comes 



^ta^ = -^lt (a2(l -/?))" + A/'' {aP)": 



(4.9) 



Since the computational grid does not extend to infin- 
ity. We adopt outgoing conditions at the outer boundary, 
that is, we assume that no radiation will enter the grid 
from large r. While this is not strictly true (there will 
be curvature back-scattering from the outgoing pulse), 
this assumption provides a reasonable computational so- 
lution. 

For the scalar field variables $ and 11, the outgoing 
conditions come from the condition on 0, namely scf) ^ 
F {s — ct) , with c — 1—2/3 being the speed of the outgoing 
waves. This means that 



• , , , s + 1 - 2/3 - 2s/3' 

$ + (1 - 2/3) + — — $ 

s 

8 + 1-2(3 + 2s/3' 



= (4.10) 



(1 - /3) (n + $) + LJ^±£0 = 0. 



(4.11) 



These equations are discretized as 

At$r + (1 - 2A*t/3r) A^c&r 



Atsf + 1 - 2^it(3f - 2Aits?A^/3f ^„ 



A*.r + i-2M./?r+2,._^rA^^^^_^ (4.12) 



and 



[(1 - /3) ($ + im + 1 ^^'^-"+^*^v tc = 0. 



(4.13) 



We can get approximate cond itions on a and K from 
their Sc hwarz schild forms (2.23) and (2.24) with the mass 
aspect ( 2.27 ) used in place of Af, since in vacuum, or 
in any region where the scalar field's self-gravitation is 
negligible, a and q should take on these forms. For 
large s, asymptotic expansion of these expressions gives: 



and 



K 



9 2m{r,t) 



+ 0{s-^). 



Thus, at large s we have, to leading order 
s (a — 1) ~ m(r, t) 

and 

s^K^ ~ m(r,t). 



(4.14) 



(4.15) 



(4.16) 



(4.17) 



We now deduce how m(r, t) behaves in the large s, weak- 
field limit. Since a ^ 1 and g — > we have 



m-47r / ($2+n^)ds. 



(4.18) 



From the condition on cj) we can see that $ ^ G{u) /s 
and n G (u) /s, where u = s — ct. Thus 



(4.19) 



mr^Sn G'' (u) du H (u) , 



that is, m is "outgoing" at large s. Therefore we get the 
following conditions for a and K^g: 



and 



s{a-l)^ H{s- ct) 
s^K\ ^H{s-ct). 



(4.20) 
(4.21) 
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These are discretized as 



V. INITIAL DATA 



A, is {a - 1)): + (1 - 2f,,p-) [s [a ~ 1))^ = 



(4.22) 



and 



At {s^K\)'l + (1 - 2^itP7) ^r" = 0. (4.23) 
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r 

FIG. 2. Sponge filter coeflicient function for A = 1.0 and 



The outgoing boundary condition reduces the am- 
plitude of reflections off the boundary, but unless the 
boundary is placed at very large r, these reflections can 
still interfere with the results of a calculation. To further 
reduce the reflections, we use a sponge filter as detailed 
in Briefly, this means that in the interior of the 

grid, we add a term to the wave equation which effec- 
tively applies the outgoing condition on a finite region 
rather than at a single radial location. For instance, we 
use the following modified evolution equation for 



$ = (/3$ + (1 - /3) n)' 

-V i>+(l-2/?)$' + 
k + 1-2(3 + 2sP' 



s + 1 - 2/3 - 2s(3' 



(4.24) 



where the coefficient function, v (r) , is given by 



v{r) = A 



(r - Ts)" {r 



,n+2 



(?l+l)(71-f 2) 



in the filtering region, Vg < r < rmax, and vanishes else- 
where. Here, A and n are parameters which can be ad- 
justed to tune the filter. Fig. ^ shows ly for A = 1.0 and 
n = 2, the values used in this work. 



We wish to examine the interactions of the black hole 
with compact ingoing pulses of scalar field. We can gen- 
erate nearly ingoing pulses using the following method. 
Let (j>{r,t) = F {u = r + t) jr. This gives = duF/r and 
(f)' — duF/r — F/r^ . For a compact pulse, we set F to a 
Guassian of the form 



F {u) = Au^ exp (- {u - cf /a'^^ 



(5.1) 



where d is an integer, and c is the radial coordinate of 
the center of the pulse. This results in scalar field initial 
data of the form 



= Arexp(- [r-cf /a"^^ 



d-l 



n = 



1 d{r-c) 
r 



2-/3 d{r-cy 



•(1-/3) 



(5.2) 



(5.3) 



(5.4) 



We s olve for /3,a, and K^e using equations (|2.6|), ( ^.10 ), 
and (|t|). Ostensibly, K\{Q,r) can be freely specified 
and the constraints can still be satisfied by appropriate 
adjustments to the other geometric variables. However, 
an arbitrarily chosen K^^[Q^r) combined with the ini- 
tial data for the scalar field given by ( |5.2| )-( ^!4| ) will not, 
in general, generate a spacetime describing the desired 
physical scenario of a scalar pulse initially infalling on a 
black hole. We therefore adopt an ansatz for i4r''r(0,r) 
which is motivated by the observation that at the initial 
time, the self-gravitation of the scalar field is generally 
relatively weak, and thus the geometric variables should 
have approximately their Schwarzschild form. Specif- 
ically, we de mand that i4r''r(0,r) satisfy the following 



equation (see (2.24)) 



-2m (r + m) 
(r [r + 2m))i 



(5.5) 



where m is the mass aspect function defined by (2.27). 
This equation, along with the equation for the mass as- 
pect (2.27) and the constraints (2.10) and ( |2.11 ), is then 
solved iteratively for any given initial scalar field con- 
figuration. We have found that initial data prepared in 
this manner does generate spacetimes of the type we seek 
even when the scalar field is significantly self-gravitating 
at the initial time. 



VI. CONVERGENCE AND STABILITY 

In order to assess the correctness and accuracy of 
our difference equations and the program which solves 
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them, we perform some tests. These include comput- 
ing convergence factors, performing a long-time vacuum 
(Schwarzschild) evolution and comparing it to the known 
solution, and performing a long-time strong-field evolu- 
tion. 

To measure convergence, we define the convergence 
factor for a grid function u by 



C 



f 



\U2h - UAh\2 
\uh - U2h\2 ' 



(6.1) 



where Ua is a solution of the finite difference equations 
on a grid with spacing a, and the £2 norm is defined in 
the usual way, that is 



\U2 



N 



N 



(6.2) 
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FIG. 3. Convergence factors as a function of time for a 
strong-field evolution with A — 2.5 x 10""^, c = 10, and a = 2. 
These values are approximately four, indicating second order 
convergence. 

The convergence factors for the scalar field and ge- 
ometric variables are shown in Fig. ^. These factors 
are computed from a strong-field evolution. They are 
approximately four throughout, indicating second order 
convergence. 

Fig. ^ shows the deviation of the evolved vacuum 
spacetime from that of Schwarzschild. We note however, 
that the spacetime being evolved is not Schwarzschild, 
rather it is a black hole in a "box" with rather ad hoc 
boundary conditions. 
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FIG. 4. The £2 norms of the differences of the computed 
vacuum geometrical variables and their Schwarzschild values 
verses time. Note that both axes are logarithmic. The outer 
boundary for this evolution was placed at r = 82M. 
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FIG. 5. The mass and the £2 norm of the 4-Ricci scalar for 
vacuum spaces with varying outer boundary position. The 
subscripts indicate the position of the outer boundary in units 
of M. The time axis is logarithmic. 

At late times, the outer boundary condition causes the 
spacetime to "drift" away from the initial configuration. 
Fig. 1^ shows plots of the mass and the Ricci scalar for 
vacuum evolutions with outer boundaries at 22M, 42M, 
and 82M. This plot clearly shows the rather large ef- 
fect the position of the outer boundary has on the late- 
time evolution. We could get better convergence of the 
difference solution to the continuum solution (i.e. with 
boundary conditions only at spatial infinity) by matching 
the interior Cauchy evolution to an exterior characteris- 
tic evolution (see [^6|jl^]) or by using an adaptive mesh 
refinement algorithm to push the the outer boundary to a 
large radius without unduly increasing the computational 
load. Nevertheless, the position of the outer boundary 
has no effect on stability. This can be seen clearly in Fig. 
^. The curves in this graph are asymptoting to fixed val- 
ues. Thus, although the computed spacetime drifts from 
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Schwarzschild, it eventually reaches a fixed configuration. 
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FIG. 6. The mass (solid line) and the norm of the 4-Ricci 
scalar (dotted line) for a vacuum evolution with outer bound- 
ary at r = 22M. Both functions are asymptoting to a fixed 
value, indicating that the spacetime is settling down to a 
static configuration. The time axis is logarithmic. 
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FIG. 7. A log-log plot of the £2 norm of the scalar field 
verses time for a strong-field evolution. The inset shows an 
expanded vertical axis for the late time evolution. The outer 
boundary is at r = 42M. 
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FIG. 8. The mass (solid line) and the I2 norm of the 4-Ricci 
scalar (dotted line) verses time for a strong-field evolution. 
The mass should be nearly constant except for a small amount 
of scalar field which radiates to infinity. The plot of the Ricci 
scalar shows the drifting of the geometry caused by the outer 
boundary conditions. Contrary to the impression given by 
this figure, the Ricci scalar is asymptoting to a fixed value. 
The time axis is logarithmic. 

Fig. shows the £2 norm of the scalar field during a 
strong field evolution to t = lOOOOM. Fig. || shows the 
mass and the Ricci scalar for the same evolution. While 
the scalar field and the mass fall off as expected, the 
plot of the Ricci scalar again shows the "drifting" of the 
geometry. 

The plots of the Ricci scalar were made by discretizing 
to first order, the analytic expression for R derived in the 
usual way by i? = where R^y is the Ricci tensor. 

Although the values of R appear large given that they 
are computed for a vacuum spacetime, they do converge 
to zero to first order in the mesh spacing as expected. 
Moreover, the values of the individual additive terms in 
the expression for R are orders of magnitude larger than 
the scalar itself, indicating that R is the size we would 
expect for first order differencing at these resolutions. 



VII. MASS SCALING 

As discussed in Sec. |^ we focus study on the evolution 
of initially ingoing "Gaussian" pulses of scalar radiation. 
In this case, the infalling field exhibits two limiting be- 
haviors, dependent on the amplitude and width of the 
pulse. These are: scattering from the existing black hole 
and collapse to form a new horizon outside the existing 
horizon. Generically, these two behaviors are separated 
by a critical value of either amplitude or width. Fig. ^ 
shows the path of the apparent horizon for various am- 
plitudes of initial data. 

The final mass of the black hole should scale as a power 
of the amplitude of the initial pulse. To fi nd ou t what 
this power should be, we can use equation (2.28). Since 
the mass is conserved, Moo is a constant. However, M is 
not constant. As the scalar field encounters the horizon, 
some mass will be transferred from the integral term to 
M. The mass of the black hole will increase by an amount 
proportional to the mass in the scalar field. By this we 
mean that after the interaction we have 



M M - 



Ank I s~ 



$2 + 

2a2 



sK' 



dr, (7.1) 



where fc is a positive constant less than 1. For a very 
narrow pulse, the entire mass of the field will go into 
the black hole, and hence, k will be very close to 1. If 
the pulse if very wide, however, k will be very close to 
0. Thus, to see how the final mass of the black hole 
scales with the amplitude of the scalar pulse, we need 



only examine the integral term in equation (7.1) 
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FIG. 9. Schematic motion of the horizon for various am- 
plitudes of the scalar field. Notice that this diagram uses the 
areal coordinate s and not the radial coordinate r. The solid 
dark vertical line which jogs right and then continues verti- 
cally represents the critical path of the horizon. The dotted 
lines are sub-critical paths and the vertical dashed lines are su- 
per-critical paths. The two thin, diagonal lines represent the 
bounds of the ingoing pulse of scalar field. A super-critical 
pulse moves inward until it crosses its gravitational radius. 
Once this happens, the apparent horizon jumps from its initial 
position to this new position where it remains. A sub-critical 
pulse moves inward until it encounters the horizon. If the field 
is very weak, the horizon is unaffected. For stronger fields, 
the horizon moves out until the pulse is entirely inside. For 
a critical pulse, the horizon moves out at the speed of light. 
Note however, that unless the energy density is a square wave, 
the horizon will not move along the straight lines as shown 
in the diagram, but will move along a curve with gradually 
increasing and then decreasing slope. 

From equations ( |).2|) -( ^!4|) we see that <i> and 11 are 
proportional to and hence to A. This means that the 
integrand is proportional to A^. This of course assumes 
the dependence of a, K^g, and /3 on (/) is much less than 
the dependence of m{r) on (p. While this seems a reason- 
able assumption, it must be checked numerically in the 
strong field regime. 
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FIG. 10. A log-log plot of the final black hole mass verses 
the amplitude of the scalar field pulse for initial data with 
d = 2, (J = 2, and c = 10. The squares represent data with 
amplitude less than the critical value, while the crosses repre- 
sent data with amplitude greater than the critical value. The 
slope of this line is 2.01 showing that the black hole mass 
depends on the square of the scalar field amplitude. 
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FIG. 11. A log-log plot of the final black hole mass verses 
the amplitude of the scalar field pulse for initial data with 
d = 4, (T = 2, and c = 10. The slope of this line is 1.99 indi- 
cating again that the black hole mass depends on the square 
of the scalar field amplitude. 



Our numerical results verify this calculation. For in- 
stance, the data in Fig. |l^ is fit by the line 

log (Af - 1) = 2.01 log A + 4.96, (7.2) 

indicating that the mass grows with the square of the 
amplitude as expected. The graph also shows there is 
no difference in behavior for sub- and super-critical data. 
That is, the final mass of the black hole exhibits the 
same dependence on the amplitude when the hole grows 
by accretion as when it forms by collapse. 
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Similarly, the data in Fig. |Tl| is fit by the line 

log(M- 1) = 1.99 logvl + 5.19, (7.3) 

indicating that the mass scaling is independent of the 
exact shape of the ingoing pulse. 

VIII. TAILS 

Fig. ^ shows <j) at constant r for runs with rmax = 
42, 82, 162. It is clear that the position of the outer 
boundary has a large effect on the late-time fall-off of the 
scalar field. Even with the sponge filter, there is enough 
reflection from the outer boundary to cause the field to 
die off more slowly than it otherwise would. 

However, with the outer boundary at rmax = 162, it 
takes at least 300M for reflections from the scattered 
pulse to travel in from the outer boundary and interfere 
with measurements at r = 30, and still longer for re- 
flections to interfere with measurements at the horizon. 
This should provide enough time to accurately measure 
the rate of fall-off of the scalar field. Fig. ^ shows the 
evolution of at r = 30 and at the horizon up until 
t = 300Af . A fit to the r = 30 curve between 200Af and 
300M shows (p falling off as t^^-^^. A fit to the horizon 
curve over the same range shows (p falling off as t^'^ '^^. 




400 



FIG. 13. A log-log plot of the absolute value of the scalar 
field at the horizon (solid line) and at r = 30M (dotted 
line) verses time for an evolution with the outer boundary 
at r = 162M. 

Linearized perturbation theory predicts these expo- 
nents should both be —3 ||l^. We note that previous 
evolutions carried out by Gundlach, Price, and Pullin 
gave exponents between —2.63 and —2.74 for at con- 
stant r HI . 



IX. RINGING 
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FIG. 12. The absolute value of the scalar field at r = 30M 
verses time for various spatial domains. Notice the differences 
in late time fall off caused by the difi'erent positions of the 
outer boundary. The vertical axis is logarithmic. 
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FIG. 14. log 101 at the horizon verses time for various pulse 
widths {A = 2.0 X 10"**). While the curves match well for 
about lOOM, there are differences in late-time behavior ex- 
hibited by the wider pulses, particularly those with a = 5 
and (T = 6. This is an outer boundary effect. As the width 
becomes larger, less and less of the initial pulse is absorbed 
by the black hole. This means there is more scalar field avail- 
able to be refiected from the outer boundary and this causes 
differences in the late-time evolution. 
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FIG. 15. log 101 at r = 30 verses time for various pulse 
widths {A = 2.0 X 10~**). 
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FIG. 16. log 101 at the horizon verses time for various pulse 
amplitudes {a = 2.0, amplitudes are xl0~*). 
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FIG. 17. log 101 at r = 30 verses time for various pulse 
amplitudes (cr = 2.0, amplitudes are xl0~^). 



Figs. |14Hl7| show the waveforms generated by the 
scattering of packets of various widths and amphtudes. 
The scalar field is measured both at the horizon and at 
r — 30M. The oscillation period in each of these fig- 
ures is approximately 53A/, and is independent of initial 
amplitude and pulse width. Gundlach et. al. measured 
an oscillation period of approximately 47Af p9| during 
collapse of Gaussian wave packets. 



X. COORDINATE EFFECTS 

Certain evolutions exhibit interesting effects which re- 
sult from the use of MMIEF co ordi nates. Recall that 
the shift comp onent, /3 satisfies ( |2.6| ). At the apparent 
horizon, ( 2.191 ) holds so we have 



2 2' 



(10.1) 



From this we can see that when no matter is crossing 
the horizon, (3 — 1/2 so the outgoing characteristic speed 
( ^.9|) is zero, as it must be since the tracked surface is 
marginally trapped. However, if / = 1, then (3 = 1 and 
the outgoing characteristic speed is —1. In this case, 
the light cone is degenerate. In fact, from ( |2.6| ) we can 
see that if / = 1, then /3 = 1 everywhere. Does / ever 
equal one? The most likely place for this to happen is 
the critical solution because that is when the "maximum" 
amount of energy is crossing the horizon for a given pulse 
shape. The values of /3 at the horizon and at the outer 
boundary are plotted in Fig. |lj for a near-critical so- 
lution, f] gets as large as .95, but never reaches 1. A 
near-critical solution from a family of pulses with d = 4 
gives a slightly higher maximum /3, but still less than 
1. It may be possible that a narrow enough pulse could 
cause / to reach 1 for an instant, but this has not been 
verified. Furthermore, it may be that such a narrow pulse 
would collapse to form a new horizon before crossing the 
existing horizon. 




10 12 14 16 



12 



FIG. 18. P at the horizon (sohd hne) and the outer bound- 
ary (dotted hne) for a near-critical solution with a = 2, d = 2, 
and c = 10. The outer boundary is at r = 42M. Outgoing 
pulses of scalar radiation will exhibit retrograde motion when- 
ever (5 > .5. 
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FIG. 19. Evolution of m' for a near-critical solution with 
a — 2,d = 2, and c = 10. Retrograde motion is apparent 
near t = 7M and t = lOM. The frames are spaced IM 
apart in time. The vertical scale changes at t = 5M so that 
the small outgoing pulse can be observed. The thin vertical 
lines passing through the frames provide a common horizontal 
reference to aid detection of the retrograde motion. 

Nevertheless, whenever /3 > .5, the outgoing character- 
istic speed is negative. Thus, outgoing pulses will appear 
to move inward when plotted in the radial coordinate, r. 
Fig. |l^ shows an evolution of m' for the critical solution 
described above. There are two periods of backwards 
motion; one at about 7M and the other at about lOM. 
These are the times when each of the "bumps" crosses 
the horizon. The retrograde motion is easier to see in 
Fig. ^ which shows contours of m' on a spacetime plot 
for the same evolution. Fig. |l] shows a fairly weak field 
evolution of m' for comparison. There is no retrograde 
motion in this case. 
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FIG. 20. Contour plot of m' for a near-critical solution with 
a = 2,d = 2, and c — 10. The small pulse moves out between 
the first and second arrows then moves back in between the 
second and third arrows. It then moves out briefly between 
the third and fourth arrows and then moves back in between 
the fourth and fifth arrows before renewing its outward mo- 
tion. 
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FIG. 21. Evolution of m' for a weak-field solution with 
A — .001, a = 2,d — 2, and c = 10. The frames are spaced 
IM apart in time. The vertical scale changes ai t — 5M the 
small outgoing pulses can be observed. There is no retrograde 
motion apparent. 

To completely remove these coordinate motion effects, 
we can abandon the shifted areal coordinate s and use 
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the usual lEF coordinates. In this case, r is areal again 
so 6 = 1 instead of 1 + //r. The evolution and constraint 



XI. NONLINEAR EFFECTS 



equat ions are trivially derived from equations (2.10)- 
( 2.16 ) with the substitutions s r and s — > 0. Note 
that the function /(<) no longer appears. 
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FIG. 22. A comparison of (f> computed in lEF coordinates 
(solid line) compared with (f) computed in MMIEF coordi- 
nates (dotted line). The frames are spaced 2M in time. The 
horizontal axis extends from to 45M in areal radius. Dif- 
ferences between the two computations are caused by small 
differences in the position of the inner boundary. 

Such a coordinate choice has the additional benefit of 
slightly simplifying the equations. However, the inner 
boundary of the grid will not be tied to the apparent 
horizon as it is in MMIEF coordinates. This causes no 
problems in practice. As matter crosses the horizon and 
the black hole grows, some of the inner grid points are 
"lost" into the black hole and we simply stop applying 
our difference equations there. Since the radial coordi- 
nate is areal, the growth of the black hole is bounded by 
the total mass of the spacetime, so there is no danger 
of having a large fraction of the grid fall into the hole. 
The previous one-dimensional black hole excising calcu- 
lations QJ^,0,^ used horizon-locked radial coordinates so 
that once the apparent horizon formed, no further grid 
points would fall into the hole. Thus, they should be 
expected to exhibit similar coordinate motion effects to 
those seen in the calculations using MMIEF coordinates. 

Fig. ^ shows a comparison of (j> evolved in IFF co- 
ordinates with (j) evolved in MMIEF coordinates for a 
strong-field case. The two match closely with differences 
caused by small differences in the positions of the inner 
boundary. 
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FIG. 23. m' at t — for a = 2, d = 2, c — 10, and various 
amplitudes, including near-critical {A*). 
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FIG. 24. m' at t = 40M for a = 2, d = 2, c = 10, and 
various amplitudes. The inset shows 1000 times fraction of 
mass scattered. For instance, the near-critical solution has 
.243 percent of its mass scattered. 

There is a sharp "bump" at the front of the outgo- 
ing pulse in Fig. This feature is absent from the 
weak-field evolution of Fig. |2l| and is certainly ampli- 
tude dependent. Fig. |23| shows a series of initial pulse 
shapes for data with various amplitudes all centered at 
r = 10. Fig. ^ shows the corresponding pulse shapes 
after scattering. 

Although this bump occurs at the front of the scattered 
pulse, it is not caused by interactions with the black hole. 
To see this more clearly, we can start the pulse further out 
and see what happens. Figs. and p6| show pulse shapes 
at t = and t — 8 ioi data with c = 20. The outgoing 
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bump develops for large amplitude data without any help 
from the black hole. A similar bump develops for initial 
data with d = A. These bumps always have the same 
characteristic shape, though their widths may vary. 
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FIG. 25. m' at t = for cr = 2, d = 2, c = 20, and various 
amplitudes (xlO"*). 
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FIG. 26. M' at i = 8 for CT = 2, d = 2, c = 20, and various 
amplitudes ( x 10"**). The small outgoing pulse appears before 
scattering with the black hole. 

It is clear from Fig. that the amplitude of this 
outgoing feature does not depend linearly on the initial 
amplitude of the pulse. For instance, from Fig. |3[ the 
height of the highest peak is about 1.5 times the height of 
the next highest peak, while their widths are the same. 
However, the amplitude of the outgoing feature in the 
first case is about 2.5 times the amplitude of the outgo- 
ing feature in the second case, while their widths are the 
same. Figures ^ and show similar nonlinear behavior 
for the pulses centered at r = 20. It is difficult to tell 



whether this feature is caused by the nonlinear interac- 
tion of an outgoing piece of the initial data with the rest 
of the pulse, or if it is caused by back-scattering from 
the effective self-potential of the ingoing pulse. Further 
study is needed. 



XII. CONCLUSIONS 

We have shown that the a null based slicing works 
well with an apparent horizon boundary condition, and 
results in a program which is stable and convergent and 
will run forever. Further, we have examined the coor- 
dinate motion effects which result from this coordinate 
system and ways to avoid them. We have also seen an 
interesting nonlinear feature in the self-gravitating scalar 
field. 

While the calculations presented in this paper were 
carried out in spherical symmetry, the coordinate system 
can be generalized to multi-dimensions. We have worked 
out an extension to the case of a Kerr black hole in 3 di- 
mensions (see |l^), the details of which will be presented 
elsewhere. 

Unfortunately, this coordinate system can not be ap- 
plied globally to spacetimes which contain more than one 
black hole since it depends on the notion of an areal 
coordinate. Nevertheless, this technique is useful for 
problems involving various matter fields and single black 
holes. For instance, it could be used to investigate the 
stability properties of the solution found by Bechmann 
and Lechtenfeld in their investigation of the scalar 
no-hair theorem (see |ll| for work in this direction). 
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APPENDIX: SPHERICALLY SYMMETRIC 
EINSTEIN-KLEIN-GORDON EQUATIONS 



In spherical symmetry (using the usual spherical co- 
ordinate names {t,r,9, 
extrinsic curvature tensor K 



the three-metric hij and the 



hij = diag 



are diagonal. We have 



{a^{t, r), r^b^it, r), r^b^ sin^ 9) (Al) 
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K'j = diag {K''r{t, r), K%{t, r), K%) (A2) The Hamiltonian constraint is 



/3* = {f3{t, r) , 0, 0) a = a{t,r) (t> = (t>{t, r) (A3) 

(A4) 

The nonzero components of the Christoffel symbols are: 



J- re 



dra rbdrirb) _ d rjrb) 

rh 

_ on rbdrjrb) . _ drjrb) 

i <A(4 = — sm U ^ i ^r<A = ; — 

rb 



(A5) 
(A6) 



r^^^ = - sin 61 cos 61 T'^^g = - cot 6i (A7) 
The two non-zero components of the Ricci tensor are 



2 „ dr{rb) 

rOr 

arb a 



a — dr\ —dr (rb) 
' a 



(A8) 



(A9) 



The evolution equations for the metric components are 
d = -aaK^'r + (a/3)' (AlO) 



b = -abK% + - (rpy . 

r 



(All) 



The evolution equations for the components of the ex- 
trinsic curvature are 



K\ = (3K''' + aK\K ■ 



1 ^a'\' 



2a 
arb 



(rb)' 



a \ a 

a-' 



(A13) 

The massless Klein- Gordon equation is formulated in 
terms of the two auxiliary fields $ and 11 which are de- 
fined by 

$ = 0' u=^(j)- /30') (A14) 
With these variables, the Klein-Gordon equation is 

$=(/3$+^ny (A15) 



2 

arb 



a J rb 



rb 



(rb)' 



(A17) 



(A18) 
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(A12) [14 



and the momentum constraint is 

rb / u ^ 
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